
CALCULATION OF AXISYMMETRIC 
SUPERSONIC FLOW PAST BLUNT BODIES 
WITH SONIC CORNERS, INCLUDING A 
PROGRAM DESCRIPTION AND LISTING 


by Jerry C. South, Jr. 

Langley Research Center 
Langley Station, Hampton, Va. 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


WASHINGTON, 




TECH LIBRARY KAFB, NM 

111 


□ 131 


75 


NASA TN D-4563 


CALCULATION OF AXISYMMETRIC SUPERSONIC FLOW PAST 

BLUNT BODIES WITH SONIC CORNERS, INCLUDING 

A PROGRAM DESCRIPTION AND LISTING 

By Jerry C. South, Jr. 

Langley Research Center 
Langley Station, Hampton, Va. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 



CALCULATION OF AXISYMMETRIC SUPERSONIC FLOW PAST 
BLUNT BODIES WITH SONIC CORNERS, INCLUDING 
A PROGRAM DESCRIPTION AND LISTING 

By Jerry C. South, Jr. 

Langley Research Center 

SUMMARY 

A program for the approximate calculation of supersonic flow of an ideal gas past 
blunt bodies with sonic corners is described. Numerical solutions are obtained for the 
system of differential equations derived from the one-strip integral (Belotserkovskii) 
method. Aerodynamic results of interest include the surface pressure and velocity dis- 
tributions, the shape and location of the detached bow shock wave, and the forebody pres- 
sure drag coefficient. Comparison of the calculations with experimental data is given, 
with particular emphasis on spherically blunted cones of large apex angle. For this con- 
figuration, a simple method is suggested for estimating angle-of-attack pressure distri- 
butions. Other shapes included are a circular disk and a spherical cap convex or concave 
to the stream. The program operation and listing in FORTRAN IV language are given in 
the appendixes. 


INTRODUCTION 

The use of aerodynamic means for deceleration of a body entering a rarified plane- 
tary atmosphere such as that of Mars brings high-drag configurations into consideration. 
In this regard, a blunt body with subsonic flow over most of its "face" is suggested; 
shapes such as a disk, spherical cap (concave or convex to the stream), and blunted cone 
with large apex angle thus assume practical importance. Under appropriate conditions, 
these particular shapes have a common feature: the transition from subsonic to super- 
sonic flow occurs in a singular manner at the sharp, convex corner which terminates the 
"face." Newtonian aerodynamic analyses do not account for the rapid pressure drop near 
the sonic corner, so other approximate methods are sought. ’Inverse" methods, where 
the shape of the bow shock wave is given and the body shape is solved for (refs. 1 to 3), 
are avoided since they can produce only analytic body shapes, whereas the shapes already 
mentioned have a singularity (corner). Among the "direct" methods (body shape given), 
one which is particularly adaptable to this problem is the method of integral relations, 
often called Belotserkovskii's method (refs. 4 to 7). In references 6 and 7, this scheme 



was specially developed for use in the first (one-strip) approximation with application to 
blunt bodies with sonic corners. 

Since there is presently a need for more information concerning the nonanalytic, 
high-drag shapes, this paper presents further results following the method of refer- 
ences 6 and 7, with primary emphasis on spherically blunted cones. The computational 
program is restricted to an ideal gas (constant specific-heat ratio) and to the singular 
sonic-point transition; this latter restriction is primarily a geometrical one. A rough, 
but general, requirement for the occurrence of the singular sonic transition at the corner 
is that the surface inclination should be everywhere greater than that at the regular sonic 
point of a full hemisphere. For example, the corner of the convex spherical cap or the 
junction of the sphere-cone must lie upstream from the usual sonic-point location for the 
full hemisphere. For any general convex shape with a corner, there is obviously a criti- 
cal combination of free-stream conditions and minimum surface angle where a smooth- 
transition sonic point will occur somewhere ahead of the corner. Difficulties are thus 
encountered when conditions are very near the critical combination, as will be discussed 
later. 

Details of the program symbols, operation, and listing in FORTRAN IV language 
are given in appendixes A to C. 


SYMBOLS 

2 

Flow variables are nondimensionalized as follows: pressure by p^V » , density 
by p^, and velocity by Voo. 

a local speed of sound 

C D drag coefficient 

D maximum body diameter 

K body surface curvature, -d0/ds 

M Mach number, V/a 

p pressure 

p m „ v stagnation pressure behind normal shock wave 

ITlcUv 

Rk base radius 


2 


nose radius 


Ri 


n 


°I 

x,r 


a 


v 

6 

$ 

X = 0 - 6 


curvilinear coordinates along and normal to body surface (fig. 1) 


value of s at sphere-cone junction 


cylindrical coordinates (fig. 1) 


velocity component in s and n direction, respectively 


total speed, \j u^ + v^ 
angle of attack 
shock angle (fig. 1) 
ratio of specific heats 

shock-layer thickness along n-coordinate (fig. 1) 
stagnation streamline isentrope constant, Pi(0)y/jpj(0^ ' 


e 

^c 

^c,det 

^c,son 

r\ 

P 


surface inclination angle (fig. 1) 
cone apex half-angle 

apex half-angle of pointed cone at shock detachment 

apex half-angle of pointed cone with sonic surface speed 

fraction of sonic speed for lower bound of surface velocity extrapolation 
bracket (eq. (7a)) 

density 

combined entropy-continuity flow variable, (p /p$)^' 



Subscripts: 


0 quantity along surface (n = 0) 

1 quantity along shock wave (n = 6) 

°° free-stream conditions 

Superscript: 

* conditions where Mq = 1 


METHOD 

Although the basic features of application of the integral method to blunt bodies with 
sonic corners have been reported before (refs. 6 and 7), there are certain minor varia- 
tions in procedure which can lead to differences in the numerical results (ref. 8). For 
the sake of completeness, then, a brief outline is given here. 

After the governing partial differential equations (in conservation form as presented 
in ref. 5) are integrated across the shock layer by approximating certain integrands as 
linear in n/6, three ordinary differential equations are obtained, as follows: 

— = (1 + K5)tan X (1) 

ds 


d/3 _ 
ds 


(1 + K6)(tan X - 5 - )p l n 1 v 1 - (2 + K5)^ p^ 2 + K8p Q u Q 2 


+ 



. r l 8 p l u l V l ) 
, r 0 / 


du Q 

ds 


^ ^ „ ,d6 6 sin 0 

7 1“1 - O u o)HJ - L 


|f ^( TlUl )df + (2 + K6)T 1 V 1 


T 0 U 0 + + K6) T jU^ 

2 ' 


6t 0 1 - M 0 - 


-1 


( 2 ) 


( 3 ) 
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A sketch of the geometry and coordinates is shown in figure 1. A particular body con- 
tour is specified by giving the surface angle and curvature, 6 and K, as functions of s. 
The functions at the shock wave are explicit functions of y, M m , j 3, and 0. The main 
dependent variables are 6, /3, and u 0 , while p Q and p Q are obtained as explicit 
functions of y, M TO , and u Q by using the isentropic law. 

The following notes are offered concerning the derivation of equations (1) to (3) : 

1. The body-oriented, curvilinear coordinates s,n were used because of their 
adaptability to a large class of smooth contours (continuous slope up to the sonic corner) 
and, in particular, to the approach to the singular sonic corner. 

2. The "entropy-continuity" formulation, so -called by Xerikos and Anderson (ref. 8), 
was used instead of the alternate "pure -continuity" formulation because the former gives 
somewhat better results than the latter in the first approximation. 

3. When integrands such as the function rur were approximated by linear inter- 
polation in n/6 between the values at the shock wave and surface, the entire function, 
including the geometric variable r, was assumed linear in n/6. 

Boundary Conditions and Solution Procedure 

At the symmetry axis, s = 0, the body surface is normal to the stream (0(0) = it/ 2), 
and the following symmetry conditions hold: 


m = tt/2 

(4) 

u 0 (0) = 0 

(5) 


At the corner of the body, the surface speed is required to reach sonic value, 

uq(s*) = a* (6) 

where a* is the critical speed, a constant which depends on y and Moo. The integra- 
tion of equations (1) to (3) starts at s = 0 and terminates at s = s*; the initial shock- 
wave standoff distance 6(0) is unknown and must be chosen so that equation (6) is satis- 
fied at the correct corner location. 

Singularity at the corner. - Inspection of equation (3) shows that the velocity gradi- 

* ^ 

ent will become infinite at s = s*, since at that point the factor 1 - Mq vanishes in the 

denominator, while the numerator does not vanish in general (with exceptions as noted 
later). In fact, as observed in reference 6, the system of equations (1) to (3) behaves 
such that the surface speed has a half-power type of singularity in a neighborhood of s* 
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a* - uq ~ const, (s* - s) 


(7) 


1/2 


while /3 and 6 are smooth (continuous s-derivative) functions at s*. An extrapola- 
tion procedure incorporating this feature is used in the numerical calculations when the 
velocity reaches a value 


77a* = uq = 0.999a* 


(7a) 


(A value of 77 which appears to give consistent results is 77 = 0.95.) 

Determination of 5(0) .- As mentioned before, the initial value of the shock-wave 
standoff distance 6(0) is unknown and must be chosen so that the singular sonic point 
occurs at the known location of the corner. Between successively improved upper and 
lower bounds, a halving-mode iteration procedure for determining the correct value of 
6(0) is easily automated; for if 6(0) is too large or too small (within certain limits), 
the sonic singularity occurs downstream or upstream from the corner, respectively. 

This behavior is not guaranteed when the desired integral curve lies near the saddle point 
of the system of equations (1) to (3), mentioned in the next paragraph. Values of 6(0) 
which are too large can yield solutions of the other family which have a behavior totally 
different from that indicated by equation (7) (the numerator of eq. (3) passes through zero 
before sonic speed is reached). However, such occurrences are easily accounted for in 
an automatic program. (See also the section entitled limitations of Program M for a dis- 
cussion of sensitivity near the saddle point.) 

It should be pointed out that the solution of the blunt-body problem by the one-strip 
integral approximation is technically easier for the singular sonic-point condition than for 
the regular sonic-point condition associated with smooth shapes without corners. In the 
latter case, the sonic point is a saddle point of equations (1) to (3) whose location is not 
known a priori (ref. 4); the desired solution is the integral curve which passes through 
the saddle point and divides two divergent families of integral curves. The inherent 
instability of the saddle-point solution makes its numerical construction difficult, and 
great accuracy in the initial value 6(0) is required as a rule (ref. 9). On the contrary, 
the singular sonic-point solutions are of one family, characterized in behavior by equa- 
tion (7); in most cases of practical aerodynamic interest, the desired solution is not too 
near the saddle point, and the relation between 6(0) and s* is nearly linear. (In fact, 
for the circular disk, the relation is exactly linear (ref. 9).) Furthermore, since the 
sonic corner lies at the maximum body radius, much of the useful information is complete 
without carrying the solution beyond the corner. 
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Remarks on Higher Approximations 

Higher approximations for the present scheme require the introduction of more 
strips between the shock wave and body surface. The regular sonic-point condition 
applies on each intermediate strip boundary, and thus the saddle points appear with a 
corresponding number of unknown initial parameters. Further technical complications 
will arise in this scheme if a portion of the subsonic flow, and hence the locus of saddle 
points, lies beyond the surface normal drawn from the body corner. In this case, polar 
coordinates might be introduced at the singular corner to complete the determination of 
the solution, but it is not known if such an approach would be successful or worthwhile. 
The difficulties associated with the higher approximations certainly justify attempts at 
obtaining as much useful information as possible from the relatively simple one -strip 
approximation. 

Before these remarks are closed, it should be mentioned that Belotserkovskii has 
reported some calculations for spherical caps (ref. 10), wherein he combined the two- 
strip approximation of the so-called scheme H of the integral method (polynomial approx- 
imation and explicit integration along the body, numerical integration from the shock 
wave towards the surface) with Vaglio-Laurin's asymptotic solution for the flow near the 
sonic corner (ref. 11). 


RESULTS AND DISCUSSION 

The program as listed in appendix C has a built-in capability to treat four basic 
axisymmetric "truncated” shapes: a disk normal to the stream, a spherical cap convex 
or concave to the stream, and a spherically blunted cone. The program can be modified 
to include other shapes if desired. (See appendix B.) In this section some typical results 
for these shapes are compared with experimental data, followed by selected results of an 
exploratory nature for the spherically blunted cone. 

Comparison With Experimental Data 

An indication of the accuracy and capability of the method can be obtained by com- 
parison with experimental data. The experimental data for the blunted and pointed cones 
were obtained in the Langley Unitary Plan wind tunnel (UPWT). The pressure data for 
both blunted and pointed cones and the shock-wave shapes for the blunted cones were 
supplied by Robert L. Stallings, Jr., Lana M. Couch, and Dorothy H. Tudor; the drag 
measurements and the shock-wave shapes for the pointed cones were supplied by 
James F. Campbell. These data are identified in the figures as UPWT data. 

Disk an d spherical caps (convex to the stream) .- Shown in figure 2 are the calcu- 
lated and measured (ref. 12) pressure distributions for a circular disk and two different 
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spherical caps (convex to the stream) for = 4.63. The calculated pressures are 
consistently too low as the sonic corner is approached, but the accuracy is adequate for 
engineering design purposes. 

Spherical cap conc ave to the stream. - In figure 3, the calculated pressures for a 
spherical cap concave to the stream are compared with those of reference 13 for 
Moo = 4.76. There is some doubt as to the plotting accuracy of the experimental data 
since it was replotted from a small-scale figure with a coarse grid (model A in fig. 6 of 
ref. 13). 

Blunt 60° cone .- The calculated and measured results for the bow shock and pres- 
sure distribution for a spherically blunted, 60° half-angle cone are shown in figure 4 for 
Moo = 2.96 and 4.63, with a bluntness ratio R n /Rb = 0.25. Here also the pressures near 
the corner are too low by, at most, 3 percent; but again, this accuracy should satisfy 
engineering requirements. The predicted shape and location of the bow shock wave is 
seen to be excellent. 

Pointed cone with d etached shock.- For y = 1.4, the bow shock is always detached 
at any value of Moo for a 60° half-angle pointed cone. Hence, when 9 C = 60° the limit 
R n /Rb — 0 presents no alteration of the important boundary conditions, and it may be 
expected that the gross features of the solutions for small values of Rn/Rb approach 
the solution for the pointed cone.* Figures 5 and 6 illustrate this feature, showing the 
calculated results for the bow shock and pressure distribution for 9 C = 60° and 
R n /Rb = 0.02 compared with experimental data for a pointed 60° cone. It can be seen 
that, at least when 9 C > 9 c ,det> the calculations for small bluntness ratios are good 
approximations to the experimental results for both the shock shape and pressures. 

On the contrary, when 8 C = 9 c ,d e t> the limit R n /Rb — 0 does involve discontinu- 
ous alteration of an important boundary condition; that is, a pointed cone (R n /Rb = o) has 
an oblique shock wave attached at the vertex, and hence the surface entropy is less than 
that corresponding to the detached, normal shock conditions which hold for any finite 
blunting (R n /Rb > Oj. Even for a pointed cone, the passage from an attached oblique 
shock to a detached normal shock at 9 C = 0c,det implies a discontinuous jump in the 
surface entropy, and hence a corresponding jump in both the vertex and sonic corner 
pressures. In spite of this, it seems that the drag of a pointed, finite cone is a continuous 
function of the cone angle while passing through 8 C = 0 c>( jet- Busemann (ref. 14) argued 

•^Sorne local nonuniform behavior is anticipated near the stagnation point as 
R n /R b — 0, since the apex of the pointed cone (with detached shock) is a singular point. 
The behavior of the flow near the apex is analogous to that of a wedge in incompressible 
potential flow, where fluid properties vary like s v , 0 < v < 1, and v depends on the 
apex angle. 
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that the drag of a sharp wedge depends continuously on the wedge angle through detach- 
ment, basing his conclusions on Guderley's potential -flow hodograph analysis (ref. 15). 
Johnston's (ref. 16) experimental results for both the wedge and cone drag at Moo = 2.45 
appear to be continuous through detachment. His investigation showed that the forebody 
pressure drag coefficient for a pointed finite cone coincides with that for the infinite cone 
when 6 C = 0 C son , departs from the infinite-cone value for 6 C > 9 C} son> anc * increases 
monotonically and continuously through detachment up to 9 C = 90°. This offers a heu- 
ristic picture, then, of the dependence of Cp on 9c for a pointed finite cone. It 
remains to consider the effects of small, but finite, bluntness. 

A rigorous discussion of the interaction between small blunting and the transonic 
singularities that arise in the pointed case for 0c, son < 0c < 0c,det ( re f- 15) is beyond 
the scope of the present report and certainly beyond the means of the coarse one-strip 
integral method. It seems natural to suppose, however, that the bluntness effects will be 
contained in an entropy layer whose thickness depends on Rn/Rb an d that in the range 
®c,son < 0c < 0c,det the drag coefficient for sufficiently small bluntness will be slightly 
smaller than that for a pointed cone. The expected difference should be caused by the 
lower pressure level near the sonic corner of the blunted cone, but the difference in drag 
should tend to zero as Rn/Rb — 0 and the bluntness -entropy effects vanish. 

The one -strip integral method does not account for vanishingly thin entropy layers; 
it "weights” the surface entropy condition more or less equally with the shock conditions. 
Thus as Rn/Rb 0, spurious entropy effects remain. Evidence of this can be seen in 
figure 7, a plot of Cj) against 9 C for a finite cone. The experimental data (Campbell, 
UPWT) were obtained for pointed cones by force measurement and subtraction of the 
estimated base drag. (The Cd values for 0 C = 50° were obtained in an earlier test, 
and were judged to be relatively higher than the others due to a different base drag cor- 
rection.) The curve corresponding to the infinite-cone solution (showing both the "strong" 
and "weak" shock branches) was obtained by calculation from the approximate algebraic 
equations derived from the one -strip integral method by setting the right-hand sides of 
equations (2) and (3) equal to zero (ref. 17). The calculated curve for the finite "pointed" 
cone solution is actually for Rn/Rb = 0.1. It was found that essentially the same curve 
is obtained for all values of Rn/Rb from 0.01 to 0.5. (See the next section and fig. 8 
concerning the effect of bluntness on drag.) It is clear that the present calculations do 
not approach the curve for the infinite cone at 6 C = 0c, son- The gap between the two 
solutions widens as Moo increases, which is in accordance with the increasing difference 
in the attached- and detached-shock entropy levels at the surface. 
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Selected Results for the Blunted Cone 

Among the shapes considered herein, the spherically blunted cone is by far the 
most versatile candidate for the aeroshell design of a planetary probe. Both the blunt- 
ness ratio and the cone angle can be varied to achieve a compromise among the require- 
ments of convective and radiative heating, drag for deceleration, aerodynamic stability, 
volume, and so forth. In this section selected results for the blunted cone are presented 
to indicate some effects of the parameters R n /Rb, &c> Y> an d M^. 

Insensitivity of drag with bluntness ratio Rn/Rb.- One unexpected result was that 

the drag remains practically constant as the bluntness ratio is changed from a nearly 
pointed cone (R n /Rb IsO.l) to a spherical cap ^R n /Rb = sec 0 C ). In figure 8 this behavior 
is shown for y = 1.4, 6 C = 60°, and = 3, 5, and 10. It is seen that for 
0 < R n /Rb < 1-5, the change in drag is not measurable, and the drag of the spherical cap 
(R n jR-^ = 2.o) is only 3 percent greater than the pointed-cone drag. On the other hand, 
the simple Newtonian approximation for the sphere-cone drag coefficient is 

C D = 2 sin20 c + (R n /R b ) 2 cos4(9 c (8) 

which predicts a 17-percent increase in passing from the pointed cone to the spherical 
cap (for 0 C = 60°). Pressure distributions for several bluntness ratios are shown in 
figure 9 for the 60° sphere-cone, and it is seen that the smaller bluntness ratios result 
in slightly higher pressures on the outboard conical skirt where the area contribution is 
the greatest. This compensating feature is not accounted for in the Newtonian theory; 
hence, the approximation of equation (8) significantly overestimates the effect of bluntness 
on drag for the sphere-cone. 

Bluntness effect on shock standoff di stance and velocity gradient.- Figures 10 and 11 
are presented to show how the bow shock standoff distance S(0)/Rb and the stagnation- 
point velocity gradient dug/ds depend on the bluntness ratio for 6 C = 60° and 
Moo = 3,5, 10, and 1000. Over the entire range of bluntness ratios from the pointed cone 
to the spherical cap, 6(0)/Rb is practically a linear function of R n /Rb. On the con- 
trary, the behavior of the velocity gradient is quite nonlinear. As R n /Rb — 0, the com- 
bination R n dug/ds approaches a finite positive limit, while R^dug/ds becomes infi- 
nite, like ^R n /R^"l. But even the qualities of the behavior near R n / R b = 0 should be 
regarded with suspicion, for, as mentioned earlier, the velocity behaves in a precise 
singular manner at the tip of a pointed cone; this is not accounted for in these calculations 
with a finite (even though small) nose radius. 

It should also be noted, according to reference 18, that the velocity gradient pre- 
dicted by the one -strip integral method for a smooth body (without a sonic corner) is 
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apparently too large. There is no reason to assume that the present results are more 
accurate. 

Effect of y on shock shape and pressure distribution .- For a given value of Moo, 
the density ratio across the bow shock increases with decreasing y, while the calculated 
shock wave moves closer to the body. Since large density ratios and thin shock layers 
are characteristic of hypersonic real-gas flows, ideal-gas calculations for low values of 
y(y — 1) are sometimes used to simulate real-gas effects (ref. 1) — at least in regard to 
the shock shape and surface pressures. The shock-wave shape and pressure distribution 
for Moo = 10, 6 C = 60°, and R n /R\ ) = 0.25 are shown for several values of y in fig- 
ures 12 and 13. Values of y lower than 1.18 cannot be considered for this combination 
of Moo, &c> an d Rn/Rb> because of the proximity to the condition where the sonic point 
of the full hemisphere would occur at the sphere-cone junction. (See section "Limitations 
of Program.") The outstanding feature of these curves is the changing character of the 
solution as conditions are approached which are appropriate to attached-shock, fully sup- 
ersonic flow for a pointed cone. For Moo = 10 and 6 C = 60°, a pointed cone will have 
an attached shock and supersonic flow at the surface for y < 1.27 (fig. 21). It is signif- 
icant that for y = 1.27 two inflection points (d/3/ds = 0) appear on the shock wave of the 
blunted cone beyond the sphere-cone junction (fig. 12). In figure 13 it can be seen that 
the pressures have a small relative minimum just ahead of the sphere-cone junction for 
y < 1.27, followed by a large recompression beyond the junction. For reference purposes, 
the locations of the shock-wave inflection points are indicated on each corresponding pres- 
sure curve. 

Traugott (ref. 19), among others, has discussed the connection between the appear- 
ance of shock-wave inflections and surface -pressure recompressions in flows over 
blunted cones. In that study, however, the flow at the surface in the recompression 
region was already supersonic, and the outgoing compression characteristics reached the 
shock wave (and caused inflections) somewhat downstream. Here the surface flow is 
subsonic, and the recompression is sensed at the shock almost immediately. For the 
lowest values of y, the solution in the nose region approaches that of a complete, smooth 
sphere, and is practically independent of the afterbody size and shape; while the solution 
beyond the junction approaches asymptotic (Rn/Rb — 0) conical pressure and shock angle 
before the expansion at the sonic corner begins to take effect. 

Effect of 6 C on shock and pressure distribution.- For fixed values of y and Moo 
the pointed sonic cone condition can also be attained by decreasing 6 C . For y = 1.4, 

Moo = 10, the sonic angle for a pointed cone is 0 c ,son = 54.6°. Figures 14 and 15 show 
that the effects of decreasing 6 C for a blunted cone (R n /Rb = 0.25) are qualitatively the 
same as those already discussed in the previous subsection, where y was decreased. 

In figure 14, the body shape was drawn for 0 C = 51°, while only the base corner positions 


11 



were indicated for the larger angles. (Note in fig. 15 that the sphere-cone junction loca- 
tions differ slightly for each value of 0 C .) Figure 16 shows how the shock standoff dis- 
tance decreases with 0 C , and the independence of the nose solution is clearly evident for 

0c < 0c, son ~ 55 °- 

Simple pressur e estimates^ for angle of attack.- It was noted that the pressure dis- 
tributions for values of 0 C larger or smaller than the nominal value, say 0 C = 60°, 
bear a strong resemblance to the experimental distributions on the windward and leeward 
generators of a 60° blunted cone at angle of attack. To investigate the possibilities of 
making some quantitative estimates, a 60° blunted cone at 5° angle of attack was con- 
sidered. Pressure data (Stallings, Couch, and Tudor, UPWT) for that configuration at 
Moo = 4.63 are shown in figure 17, together with the present results for 6 C = 65° 
and 55°, approximating the windward and leeward pressure distributions, respectively. 

It appears that the levels of the predictions are roughly in accord, but some coordinate 
shifting and stretching might help to account for (1) the displaced stagnation point and 
(2) the altered distance to the sonic corner. One quick and obvious scheme is to shift 
the stagnation point to the most forward point of the nose and cause the sonic corners to 
occur at the correct location by a linear transformation. That is, the simulated wind- 
ward (0c = 65°) and leeward (0 C = 55°) pressure calculations are relocated by correcting 
the corresponding s/D values by a term which depends linearly on s. For this 
example (60° blunted cone at 5° angle of attack), 


and 


s_ 

wind 






0.011 - 0.012 



65 




0.001 


0.018 



55 


where A is the displacement required to shift the stagnation point to the most forward 
nose point (which lies on the spherical portion for a = 30°), 


A 

D 


aR n 

D 
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The results with the shifted and stretched coordinates are replotted in figure 18, together 
with similar calculations and data for a = 10°. Figure 19 shows the comparison for 
Moo = 2.96. The agreement for the leeward pressure distribution is surprisingly 
accurate, while the windward results agree within 5 percent. This approximation may be 
particularly useful if it is coupled with some empirical fit for the peripheral pressure 
distribution such as the second-degree trigonometric polynomial suggested by Kaattari 
(ref. 20). Hence there exists the capability of obtaining, quite cheaply, reasonable esti- 
mates for the entire forebody pressure distribution. 

LIMITATIONS OF PROGRAM 
Smooth Sphere Sonic Point at Junction 

The main difficulty in this program is encountered for the sphere-cone shape when 
conditions are approached such that the sonic point for a complete, smooth sphere would 
occur at the sphere -cone junction. ^ The sonic point for the smooth sphere is a saddle 
point of the differential equations (1) to (3) (ref. 4), and the diverging integral curves near 
the saddle point are extremely sensitive to the value of 6(0). For combinations of y, 
Moo, and 0 C near this condition, the success or failure of convergence depends on three 
factors: (1) the number of figures carried in the calculations (i.e., the "word" length of 

the computer), (2) the bluntness ratio R n /Rb, and (3) the prescribed accuracy criterion e 

for placing the sonic singularity ^JrQ* y /R] C) - l| = ej. An example follows. 

Figure 20 is an illustration of the sensitivity of the solution on 6(0) near the criti- 
cal condition. For y = 1.4 and Moo = 10 the sonic point for a smooth sphere occurs at 
a location where the surface angle 0, and hence the critical cone angle 6 C is about 
49.1°. But even for 9 C = 52°, eight correct figures in 6(0)/R n were required to place 
the sonic singularity slightly beyond TQ*/R n = 4.0. In other words, if only eight deci- 
mal figures were carried in the calculations, bluntness ratios less than about 0.25 could 
not be obtained, and those somewhat greater than 0.25 would be possible only for a fairly 
lax convergence criterion (say e ~ 0.05). For 6 C - 51°, R n /R^ = 0.25, forty-five 


^In an exact analysis, the chief concern would be the point of origin of the limiting 
characteristic (ref. 1, pp. 202-203), downstream from which changes in the body shape 
cannot affect the nose solution. In the one -strip approximation, however, it is the surface 
sonic point that is crucial. 
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halving-mode iterations (producing about fourteen correct figures^ in S(0)/R n ) were 
required for an accuracy criterion e = 10~^. 

In a study of blunted cones of "small" angle (i.e., with sonic point always on the 
spherical nose), Traugott (ref. 19) described a range of cone angles in which the integral 
method should fail, for a given y and M^. The upper limit is the critical angle just 
described; that is, the cone angle 6 C for which the sonic point of the complete smooth 
sphere would occur at the sphere-cone junction. The lower limit is the cone angle for 
which the asymptotic, pointed-cone pressure is equal to the sonic pressure obtained by 
isentropic expansion from the normal shock conditions. The cone angles in between con- 
stitute what Traugott called the "second sonic point" region; wherein the flow along the 
surface will become supersonic on the spherical nose, but must become subsonic again 
if asymptotic conical pressure is to be reached downstream. While Traugott found the 
method to fail upon approaching this range of angles from below, the present program 
fails upon approaching from above. 

Unfortunately, there are not available any extensive tabulations of the sonic-point 
location on a sphere for various combinations of y and In fact, it is that informa- 

tion for the one-strip integral method that is needed here, for the purpose of avoiding the 
critical condition when using this program. A more practical limit, which can be charted 
profusely, is the sonic condition for the pointed cone, 0c,son(^>^°°)- K was found that 
this condition lies fairly close to the "smooth sphere sonic-point" condition and will serve 
as a guide in warning the program user of expected sensitivity and possible nonconver- 
gence. Figure 21 is a chart of the pointed sonic cone condition; the calculations were 
performed by using the algebraic equations which approximate conical flow by setting the 
right-hand sides of equations (2) and (3) equal to zero (ref. 17). 

Convex or Concave Spherical Caps 

The same difficulties are to be expected for the spherical cap (convex to the stream) 
when the corner is near the smooth sphere sonic point. For the concave spherical cap, 
no extensive investigation was undertaken to determine the limits of applicability of the 
present program. An obvious geometrical condition to be avoided for this shape is the 
intersection of the surface normals within the shock layer — that is, since 

3 

A reliable estimate for the number of halving-mode iterations Nh m required to 
produce decimal figures in S(0)/R n is 

N hm “ N d/ lo gl0 2 * 3 - 3N d 
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K = -d0/ds = -1 for the concave spherical cap, the condition K6 £ 1 might occur. It 
is possible that this could happen for low values of M 1JO and/or values of R n /R b 
approaching 1.0 (from above). 


Moo Too Small 


As Moo decreases, the shock wave moves farther away from the body and the 
sonic point on the shock wave moves up, away from the axis. For a given body shape, 
a value of M^ > 1 will be reached where a considerable portion of the subsonic flow 
region lies beyond the normal drawn from the sonic corner of the body. It would seem 
that such solutions would be either nonexistent or meaningless in the one-strip integral 
approximation, since the calculation ends when the sonic corner is reached, and hence a 
part of the subsonic flow would be ignored. Such a situation is, of course, both physically 
and mathematically objectionable. 


What appears to happen in the one -strip approximation, at least in all cases studied 
so far, is the following: the shock-wave sonic point always occurs ahead of the surface 
normal drawn from the sonic corner of the body, which implies proper closure of the 
subsonic region. But since the shock wave has less curvature near the axis as Moo 
decreases, it must turn more rapidly as the sonic corner is approached to achieve the 
sonic shock angle; the factor 9 (Pi u i v i)/ 9 / 3 in the denominator of d/3/ds becomes 
small near the corner, resulting in the rapid increase in magnitude of d/3/ds. Ultimately 
a value of M^ is reached where 8 (Pi u i v i) / / a / 3 passes through zero at the corner, and 
hence the shock has infinite curvature there. The disk provides a good example of this 


difficulty. In table I are given values of 8 (Pi u i v i)/ 8 /3 as M oo decreases from hyper- 
sonic values; it is clear that two -fold smooth (continuous curvature) shock-wave solutions 
do not exist below M^ ~ 2.3. 


TABLE I.- APPROACH TO SHOCK CURVATURE SINGULARITY FOR A DISK 


M TO 

6(0)/R b 

8 (pl u l v l)*/ 8 /3 

10.0 

0.470 

2.004 

6.0 

.494 

1.812 

4.0 

.541 

1.464 

3.0 

.607 

.997 

2.4 

.703 

.284 

2.3 

.732 

.024 
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CONCLUDING REMARKS 


The present program has the capability for calculating good approximations to the 
shock-wave shape and surface pressure distribution for blunt bodies with sonic corners. 
Typical results for the disk, concave or convex spherical caps, and sphere-cones agree 
well with experimental data. 

The flow past a pointed cone was approximated by calculations for a sphere-cone 
with small bluntness ratio. It was found that results were generally good for cone angles 
larger than the detachment angle; but results were not satisfactory as the cone angle was 
decreased to and below the detachment angle. The reason for this appears to be the 
inability of the method to account for a vanishingly thin entropy layer. 

For the blunted cone, inflection points occur in the shock wave when the sonic con- 
dition for a pointed cone is reached. This condition can be reached, for example, by 
decreasing the cone angle. Further decreases in cone angle cause the nose and afterbody 
solutions to become nearly independent of each other. 

It was found that reasonable estimates could be made for angle-of-attack pressure 
distributions in the symmetry plane of a blunt cone. The windward and leeward distri- 
butions were simulated by adding and subtracting the angle of attack from the cone angle, 
respectively. 

A critical condition for the method is evident in the application to sphere-cones and 
spherical caps. It occurs when the combination of parameters (cone angle or surface 
inclination angle at the sonic corner, Mach number, and specific heat ratio) is such that 
the natural sonic-point location for the complete sphere lies at the sphere-cone junction 
or the spherical cap corner. Depending on the bluntness ratio, extreme sensitivity and 
nonconvergence is encountered somewhat before this condition is reached. 

For a given body shape, there is a minimum Mach number below which no two-fold 
smooth (continuous curvature) shock solutions exist. The accuracy of the solutions near 
this minimum value is questionable, particularly near the sonic corner. For the disk, 
the minimum Mach number was about 2.3 (in air). 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., February 13, 1968, 

129-01-03-06-23. 
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APPENDIX A 


PROGRAM SYMBOLS 


The FORTRAN symbols appearing in the input, headings, and output of this program 
(appendix C) are given in the left-hand column, and their meanings defined in terms of 
standard notation or the symbols defined in the section "Symbols" are given in the right- 


hand column. 

These program symbols are as follows: 

Input: 


LC 

body shape trigger, 1, 2, 3, or 5 for disk, convex spherical cap, blunted 
cone, or concave spherical cap, respectively 

MM 

loop counter, number of cases (combinations of y, Moo, 6 C , and R) 
for a given body shape (LC) 

GAMMA 

y 

AMINF 

Moo 

THETAC 

0 C , cone half-angle in degrees for blunt cone (LC = 3); dummy input 
value for other shapes (other LC) 

R 

sonic-point radius parameter, 1.0 for disk (LC = 1) and rQ*y/R n for 
other shapes (other LC) 

Case heading (other than symbols defined for input): 

DELO 

6(0), initial standoff distance in appropriate length scale (as explained 
in appendix B) 

DS 

initial integration step size 

ETA 

ri (see eq. (7a)) 

SI 

value of s at sphere-cone junction (LC = 3); dummy value for other 
shapes 
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NiS number of integration steps in final run 

NiR number of integration runs to find correct 6(0) 

CD forebody pressure drag coefficient (referred to base area 7 itq *2 for 

all shapes) 

Output (the output is printed across in a block of three rows for each integration step): 


s 

S 

BE TAD 

j3, degrees 

THETAD 

6, degrees 

XO 

x 0 

RO 

r 0 

UO 

u 0 

PO 

Po 

RHOO 

p 0 

TO 

temperature ratio, Tq/T, 

DELTA 

6 

XI 

X 1 

R1 

r l 

VI 

V 1 

U1 

U 1 

PI 

Pi 
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RHOl 

Pi 

DS 

integration step size 

DBDS 

djS/ds 

DUODS 

duo/ds 

POBAR 

Po/Pmax 

SBAR 

s/s* 


Main program (some of the more impor tant FORTRAN symbols used in the main 
program): 


DSC starting step size for each integration run 

LN overall print trigger, 1 for no print and 0 for print (switched to zero 

when correct 6(0) is found) 

N step print trigger, 0 for no print and 1 for print (only print when stable 

step is completed) 

L body shape trigger, equal to LC but switched to 4 for rescaling when 

LC = 3 (see description of body shape subroutine in appendix B) 

LL trigger for double-valued curvature at sphere-cone junction, 1 for 

sphere curvature and 2 for zero curvature 

LLL step size trigger, 1 until step size is halved to achieve velocity bracket 

and 2 thereafter (see program statement 54) 

LLLL halving mode trigger, 1 until upper and lower bounds for 6(0) are 

found and 2 thereafter 

Body shape subroutine : 

THETA 6, radians 

AK K, surface curvature 
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PROGRAM OPERATION 

See main text section entitled "Limitations of Program" before operating this 
program. 


General Remarks 

The FORTRAN IV program listed in appendix C was originally developed for the 
IBM 7094 electronic data processing system, and subsequently it was modified for use in 
the Control Data 6600 computer system. Other than using control cards which are 
appropriate to an individual system, the changes necessary to operate the program on the 
IBM 7094 system are as follows: 

(1) Remove the two instructions in statement 1, 

IF (EOF, 5) 9999,1000 
9999 STOP 

(2) In statement 2, change the thirteenth instruction to: 

IF (NIR.LT.31) GO TO 2000 

The reason for change (2) is the shorter single -precision word length of the IBM system 
(eight decimal digits). More than 27 halving iterations for 6(0) will not change the 
eight decimal digits carried in the calculations. About four runs are allotted for finding 
upper and lower bounds for 6(0). 

A large number of comment cards were included to highlight the various sections 
of the program. A fourth-order Runge-Kutta integration scheme was built into the pro- 
gram so that a "library" routine would be unnecessary. Also, a rule-of-thumb step size 
variation scheme is used which is very simple yet adequate for these calculations. It 
tests on the shock-angle derivative dj3/ds and seeks a step size such that the "pre- 
dicted” and "corrected" values for d/3/ds at the interval midpoint are within a few per- 
cent of each other. Numerical stability and accuracy are not a problem in these 
calculations. 

Initial estimates for DE LO.- For each of the four shapes included, the initial esti- 
mate for 6(0) is automatically computed, as follows: 
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LC 

Shape 

DELO 

Source 

1 

Disk 

1.03^ Pl (0) - 1 

Ref / 21 

2 

Convex sphere cap 

.gqi/Ip^o) - g 

Ref. 21 

3 

Sphere -cone 

.667R/[ Pl (0) - g 


5 

Concave sphere cap 

1.03r/^ Pi (0) - 1 



Accuracy criterion . - The accuracy criterion for locating the sonic corner is listed 

as 10“^; that is, |ROST/R - l| = 10”^ (see the fourth instruction from the end of state- 
ment 60). This can be relaxed, if necessary, to achieve convergence closer to the 
"critical condition" described in the main text in the section "Limitations of Program." 

Modifications for other body shapes .- The program can be modified to include other 
shapes in the body shape subroutine (BSR) by adding more statements, triggers, and, if 
there are more than two segments to the shape, more junction locations SR, Sin, etc. 

Care must be taken to insure that an integration step coincides with each junction loca- 
tion, and that the curvature of the segment upstream or downstream from the junction is 
properly used. Some study of the logic for the sphere-cone shape (LC = 3) should help 
in this respect; the crucial locations in the main program are the last "IF(. . .)" instruc- 
tions in statements 30 and 52. 

Some reasonable initial value for DELO must be given for any new body shapes, or 
otherwise the initial DELO will be estimated for a sphere. 

Input 

The first four cards (ER, BLK, BCK, ERR) are always necessary to the program, 
and they provide for specific error messages or debug information. 

NNN is the debug trigger. If NNN = 0, debug information is printed out; a dummy 
fixed-point number different from 0 is read in if no debug information is needed. 

LC is the body shape trigger; MM is the number of cases to be run for that value of 
LC. For each case, one card is required, with values for GAMMA (y), AMINF (Moo), 
THETAC (6> c ), and R. 

For LC = 3 (the sphere-cone) THETAC is the cone half-angle in degrees; for 
other LC, a dummy value must be input. 

R is the sonic corner location parameter. For LC = 1 (disk), R = 1. For the 
three other shapes, R = rQ*/R n . Note that for the sphere-cone, R is the reciprocal of 

the bluntness ratio (Rn/Rb) ■ 


21 



APPENDIX B 


Case Heading and Output 

See appendix A for a description of the symbols in the case heading and output. For 
each integration step, a three-line block of 21 values is printed. 

Length scale.- The appropriate length scale is as follows for the four body shapes: 

Disk Rj, = 1 

Concave or convex sphere cap R n = 1 

Sphere-cone Rb = 1 


Sample Cases 

At the end of the listing in appendix C, some sample input cards are listed. After 
the four error- and debug-message cards, the NNN value is 1 (for no debug). Then all four 
body shapes are run: one case each for the disk (LC = 1) and the concave spherical cap 
(LC = 5) and two cases each for the convex spherical cap and the sphere-cone. The total 
central processing unit time for all six cases on the Control Data 6600 computer system 
was about 18 seconds. The correct value of DELO (6(0) to the appropriate length scale 
as already described) is listed for each case in order in the following table (15 decimal 
figures carried in calculations): 


Shape 

LC 

MM 

GAMMA 

AMINF 

THETAC 

R 

DELO 

Circular disk 

1 

1 

1.4 

4.63 

ai.O 

1.0 

b 0. 519443847 

Spherical cap, convex to stream 

o 

9 

1.4 

4.63 

a 1.0 

0.25882 

c 0. 107979173 

Ci 

Ci 

1.4 

4.63 

a 1.0 

0.5 

CO. 156923693 

Spherically blunted cone 

Q 

9 

1.4 

4.63 

60.0 

4.0 

b 0. 0921382953 

O 

Ci 

1.4 

4.63 

60.0 

50.0 

b 0. 0556326561 

Spherical cap, concave to stream 

5 

1 

1.4 

4.76 

a 1.0 

0.25882 

c 0. 160711018 


a Dummy input. 

b S(0)/R b . 

C 6(0)/R n . 
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PROGRAM LISTING 
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c truncated blunt body 

COMMON ER ( 1 2 ) ♦ BLK ( 12 ) *BOC( 12 ) * ERR (12) * THET A * THET AD ♦ AK »RBET AD • ALA MB 
1 * CON 1 • CON2 ♦ CON3 ♦ CON4 * CONS* CON6 « C0N7 % C90 * C91 *NNN*U1 *Vi *P1 • GAMMA *RHO 
21 tSTl *Z1 #H1 *R1 «X1 ♦CPHl % RHOO « PO • AOSQ »STO*PHl * DDDS • DBDS « DUODS • DRODS t 
3DX0DSt TO ♦ AOST * BET A ♦ FAC ♦ DENOM « OENOM 1 • DFDS • CD « R *PU1PB t POBAR ♦ SMAX ♦ 
4SBAR f LN»DELTAL»DEL0f DELTAUf NlRtNl S*LtPZlP3 

800 FORMAT(25H DA V 1 S-DLD-A 1 054C-3- 2 5-66 ) 

801 FORMAT ( 1 2H JERRY SOUTH) 

802 FORMAT ( 21 H TRUNCATED BLUNT BODY//) 

805 FORMAT ( 31 H CIRCULAR DISK NORMAL TO STREAM//) 

806 FORMAT ( 31 H SPHERICAL CAP CONVEX TO STREAM//) 

807 FORMAT ( 43H SPHER I CALL Y-BLUNTED CONE*CONE HALF-ANGLE =tE16.8//> 

808 FORMAT ( 32H SPHERICAL CAP CONCAVE TO STREAM//) 

809 FORMAT (25H ITERATION COUNT EXCEEDED/) 

90 FORMAT ( 12A6) 

91 FORMAT ( 1015) 

92 FORMAT (5E14*8) 

93 FORMAT ( 1 X6HGAMM A=E 1 6 • 8 ♦ 2X6HAM I NF=E 1 6* 8 ♦ 2X5HDEL0=E 1 6 • 8 « 2X3HDS=E 1 6 • 8 
1 ) 

94 FORMAT <1X1 A6/ ( 8E I 6 • 8 ) ) 

95 FORMAT ( 1X12A6// ) 

96 FORMAT (15H DELO TOO LARGE) 

97 FORMAT ( 8X 1 HS 1 4X5HBETAD 1 0X6HTHET ADI 2X2HX014X2HR01 4X2HU014X2HP0 13X4 H 
1 RHOO/7X2HTO 1 4X5HDELTA 1 2X2HX 1 14X2HR1 14X2HV1 1 4X2HU1 14X2HP1 13X4HRH01/ 
27X2HDS 1 4X4HDBDS 1 2X5HDU0DS 1 1 X5HPOB AR 1 1 X4HSBAR// ) 

98 FORMAT (/) 

99 FORMAT ( 1 X8E 1 6 * 8/ 1 X8E 1 6 • 8/ 1 X5E l 6 *8/ ) 

100 FORMAT (//) 

1 02 FORMA T( 1 X4HETA = E 1 6 • 8 1 2X2HR = E 1 6 • 8 * 2X3HS I = E 1 6 « 8 • 2X4HN I S= I 5 * 2X4HN I R= I 
15f 2X3HCD=E16.8//) 

WR I TE ( 6 * 800 ) 

WR I TE ( 6 • 80 1 ) 

WR 1 TE ( 6 ♦ 802 ) 

DEG=57 • 295780 
PI=3. 1415927 
P I 2=P I /2 • 

READ { 5 ♦ 90 ) ER *BLK ♦ BCK * ERR 
READ ( 5 ♦ 91 ) NNN 
1 READ {5*91 )LC*MM 
I F (EOF * 5 ) 9999 * 1000 
9999 STOP 

1000 DO 71 MMM= 1 ♦ MM 

READ (5* 92 ) GAMMA ♦ AM I NF * THET AC • R 
C COMPUTATION OF CONSTANTS AND STARTING PARAMETERS 
ET A= • 95 
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5 I = • 1 E+ 06 

IF (LC*EQ*3 )S 1= ( 90 • — THET AC ) /DEG 
DSC= .05 

IF(LC»EQ»3) DSC = S I / 1 O • 

C0N3 = GA MM A — 1 • 

CON4=GAMMA+l • 

CONI =4#/C0N4 
CON2 = AM I NF ** 2 
CON5=CON3/ (2#*GAMMA ) 

C0N6=1 •/(GAMMA*C0N2) 

C0N7=1 ./C0N3 
L = LC 

DENS I *2 • * ( C0N2— 1 • ) / < 2 • + CON3*CON2 ) 

DELO= • 667 /DENS I 
IF (L.EQ.3 >DELO=R*DELO 

IF(L.ECUl#OR*L.EQ*5 )DEL0=1 .03*R/SQRT (DENS I ) 

DELT AL=DEL0 
N I R = 0 
LN= 1 
LLLL=1 

C INITIAL VALUES 
2 LL = 1 
LLL= 1 
N = 0 

DS=DSC 
S = 0*0 

CALL 8 SR (L »LL * S I « S ) 

I F ( LN»EQ #0 • AND# L.EQt 1 ) WRITE ( 6 % 805 ) 

I F (LN»EQ* 0 • AND* L»EQ • 2 ) WR I TE ( 6 * 806 ) 

IF (LN»EQ*0 • AND • LC • EQ • 3 ) WR I TE ( 6 f 807 ) THET AC 
I F ( LN» EQ • 0 • AND • L* EQ • 5 )WRITE (6* 808 ) 

IF ( LN« EQ • 0 > WRI TE (6, 93 ) GAMMA , AM I NF * DELO * DS 
IF ( LN* EQ • 0 )WRITE ( 6 « 102 ) ET A » R « S I , NI S*NIR t CD 
IF ( N I R • LT *50 )G0 TO 2000 
IF (LN*EQ*0 > WRI TE(6« 809 ) 

I F ( L • EQ »2 ) THET ST = SST*DEG 

IF (LN.EQ.O ) WRI TE (6* 94 )ER( 1 ) * DELTAL • DELO * DELT AUi SST « ROST * THET ST 
IF (LN.EQ.O >STOP 
L N = 0 
GO TO 2 
2000 NI S=0 

N I R=N I R+ 1 

IF (LN.EQ.O ) WRI TE ( 6* 97 ) 

U0=0*0 
BETA = P I 2 
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BETAD=BETA*DEG 
DELT A = DELO 
RO=0*0 
X0=0*0 
F = 0 • 0 

200 CALL SHOCK ( S *U0 ♦ BET A « DELT A * RO * XO * F ) 

D I SC = — V 1 

IF (NNN*EQ*0 ) WR I TE (6*94 ) BLK C 1 ) « ALAMB ♦ 8 * A * P 1 •RH01 ♦ V 1 * D I SC « DELTA < BETA 
1 • UO * RO ♦ XO * THETA * R 1 *X1 *Ul 
IF (DISC*LT*0*0.0R.RH01 *LE* *00001 )G0 TO 20 
AO ST = SQRT (DISC ) 

TEST=ETA*AOST 
INITIAL DERIVATIVES 

I F ( NNN • EQ • 0 ) WRI TE(6*94 ) BLK (2 ) ♦ AOST * TEST * CPH I *H1 * FAC • RHOO * PO ♦ Pu 1 P 

IB* AOSCUTO 
GO TO 24 

20 WR ITE (6 *94 )ER< 1 ) * DI SC *RH01 
GO TO 71 

24 DEN0M=DELTA*H1 *PU1 PB 

IF CABS (DENOM ) *LE. *00001 )G0 TO 26 

DBDS = ( PO— P 1 — V 1 *H 1 ) /DENOM 

DEN0M1 =DELTA*RH00*V1 

IF (ABSCDEN0M1 ) *LE* *00001 )G0 TO 27 

DUODS=- ( 1 • +AK*DELT A )* (PO-P1 )/DENOMl 

DDDS=0 • 0 

DRODS= 1 .0 

DXODS=0 • 0 

DFDS=0 • 0 

CD = 2* * ( PO— C0N6 > 

HUO=UO 

IF (LN*EQ*0 ) WRITE (6* 99 ) S * BETAD • THET AD * XO * RO * UO * PO * RHOO * TO * DELTA ♦ X 1 * 
1R1 * VI * U 1 * P 1 * RHO 1 * DS*DBDS*DUODS*POBAR, SBAR 

25 1 F ( NNN • EQ • 0 ) WRI TE (6 *94 )BL< (3) « DENOM « DBDS * DE NOM 1 * DUODS * DDDS * DRODS * D 
1 XODS 

GO TO 31 

26 WR I TE C 6 * 94 ) ER ( 4 ) * DENOM 
GO TO 71 

27 WR I TE < 6 * 94 > ER ( 5 ) * DENOM * DENOM 1 
GO TO 71 

RUNGE— KUTT A INTEGRATION (FOURTH ORDER ROUTINE) WITH VARIABLE STEP SIZE 
30 IF CL*EQ.2.AN0*RS*GE* 1 *5707 )GO TO 71 
CALL BSR<L*LL«SI * RS ) 

CALL SHOCK < RS ♦ RUO ♦ RBET A « RDELT A « RRO « RXO « RF ) 

IF(L*GT.l . AND*PZ1 PB*LT. • 1 )GO TO 530 

IF ( LN* EG • 0 • AND • N* EQ • 1 ) WRITE (6*99 )RS * RBETAD * THE TAD *RX0 * RRO ♦ RUO , PO ♦ 
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1 RHOO * TO .RDELT A * X 1 »Rl,VHUliPl ,RH01 »DSi DBDS * DUODS ♦ POBAR • SBAR 
N = 0 

IF { S 1 N ( RBET A )**2 • LE • .00001 . OR . ABS ( COS ( AlAMB ) ) . LE . .00001 •OR.FAC.LT 
1 0 • 0 • OR . ABS (RHOO ) • LE • • 0000 1 • OR. ABS ( RRO ) • LE • .00001 • OR • ABS ( DENOM } .LE 
2.00001 .OR. ABS (DENOM1 ) .LE. .00001 )G0 TO 71 
IF (LL.EQ.2 )GO TO 31 

CHECK ON LOCATION OF SPHERE— CONE JUNCTION 
IF (S+DS.LE.SI )GO TO 31 
DS = S I — S 

31 Al=DS*DUODS 
B1 =DS*DBDS 
C1=DS*DDDS 
D1 =DS*DRODS 
El =DS*DXODS 
FI =DS*DFDS 

I F (NNN.EQ.O > WR I TE < 6 .94 ) BLK C4).Al.Bl.ClfDl«El«Fl 
RS=S+0S/2. 

RUO=UO+Al /2. 

IF (RUO.GT. *999*AOST ) GO TO 54 
RBETA=BETA+B1 /2 • 

RDELTA=OELTA+Cl /2. 

RRO=RO+Dl /2 • 

RXO=XO+El /2 . 

RF=F+F 1 / 2 . • 

CALL BSR (L.LL.S I .RS ) 

CALL SHOCK (RS.RUO.RBETA . RDELT A . RRO * RXO iRF ) 

IF < S I N ( RBET A ) **2 • LE . .00001 .OR. ABS ( COS (AlAMB ) ) .LE. .00001 . OR .FAC • LT 
1 0.0. OR. ABS (RHOO ) .LE. • 0000 1 .OR. ABS (RRO) . LE • .00001 .OR. ABS (DENOM ) .LE 
2. 0000 1 .OR. ABS ( DENOM1 ) .LE. .00001 )G0 TO 71 
A2=DS*DUODS 
B2=DS*DBDS 
C2=DS*DDDS 
D2=DS*DRODS 
E2=DS*0XQDS 
F2=DS*DFDS 

IF (NNN.EQ.O ) WR I TE (6. 94 ) BLK ( 5 ) , A2 % B2 . C2 . D2 • E2 . F2 
RUO=UO+A2/2. 

IF (RUO.GT. .999*AOST )G0 TO 54 
RBET A =BET A +82/ 2 . 

RDELT A = DELT A+C2/2 . 

RRO=RO+D2/2. 

RXO=XO+E2/2. 

RF=F+F2/2. 

CALL BSR ( L * LL • S I • RS ) 

CALL SHOCK ( RS . RUO . RBET A . RDELT A . RRO . RXO . RF ) 
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IF (SIN (RBETA ) **2 • LE • • 0000 I * OR • ABS ( COS ( ALAMB ) ) *LE. .00001 *0R*FAC*LT* 
1 0*0*OR*ABS (RHOO > *LE* *00001 *OR. ABS (RRO ) *LE. *00001 • OR • ABS ( DENOM ) * LE * 
2 • 0000 1 • OR • ABS ( DENOM 1 )*LE**00001 )GO TO 71 
A3=DS*DU0DS 
B3-DS*DBDS 

IF ( ABS (B2 ) *LT* • 1 E-05 . OR .LLL • GT * 1 )GO TO 32 
C INTEGRATION ACCURACY STEP-SIZE TEST 
DSTES T = ABS ( (B2-B3 ) /B2 ) 

IF (DSTEST*LE. .05 )GO TO 32 
DS= • 5*DS 

IF (DS*GT* • 1 E-05 )G0 TO 55 

WRITE (6 *94 )ER(5 ) * DS * DSTEST * B1 . B2 * B3 * DEL T AL * DELO * DELT AU 
WR I TE ( 6 * 9 1 )NIR*NIS 
GO TO 71 
32 C3=DS*DDDS 
D3=DS*DR0DS 
E3=DS*DXODS 
F3=DS*DFDS 

IF (NNN.EQ.O ) WR I TE (6 *94 )BLK (6 > , A 3 * B3 * C 3 * D3 * E3 * F3 

RS=S+DS 

RU0=U0+A3 

I F ( RUO • GT • • 999* AO ST ) GO TO 54 

RBETA=BETA+83 

RDELT A = DELT A+C3 

RR0=R0+D3 

RX0=X0+E3 

RF = F+F 3 

CALL BSR ( L ♦ LL ♦ S I ♦ RS ) 

CALL SHOCK ( RS * RUO « RBETA * RDELTA • RRO * RXO * RF ) 

IF (SIN(RBETA )**2.LE. .00001 • OR. ABS ( COS (ALAMB ) ) *LE* *00001 *OR*FAC*LT * 

1 0.0 .OR • ABS (RHOO ) . LE* • 0000 1 • OR. ABS (RRO ).LE*. 00001 .OR* ABS ( DENOM ) *LE . 
2 • 00001 • OR * ABS ( DENOM 1 ) • LE* *00001 ) GO TO 71 
A4=DS*DU0DS 
B4=DS*DBDS 
C4=DS*DDDS 
D4=DS*DRODS 
E4=DS*DXODS 
F4=DS*DFDS 

IF (NNN.EQ.O ) WR I TE ( 6 * 94 ) BLK ( 7 ) , A4 * B4 * C4 * D4 * E4 • F4 
U0 = U0+1 •/€>•* ( A 1 +2 .*A2+2.*A3 + A4 ) 

C TEST FOR VELOCITY DECREASE 

IF(ABS(UO).LE..00001 )G0 TO 530 
I F ( UO ) 530 * 530 * 5 1 

C TEST FOR VELOCITY BRACKET NEAR SONIC POINT 
51 IF (ABS(UO-TEST ) .LE. .00001 )G0 TO 53 
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IF (UO-TEST >52*53*53 

52 IF <S.EQ*0#0 )DSC = DS 
S = S + DS 

HUO=UO 

BET A = BETA+ 1 . /6 • * ( B1 +2 • *B2+2 • *B3+B4 ) 

DELTA = DEL T A+ 1 . /6* * ( C 1 +2 * *C2+2 » *C3+C4 ) 

RO=RO+l */6»* ( D 1 +2 • *D2+2 • *D3+D4 ) 

IF (RO.GT.R*AND*L*GT. 1 » ANDiLN*EQt I > GO TO 530 
XO=XO+l •/&•* (El +2**E2 + 2 .*E3+E4 ) 

F=F+ 1 • /6 • * ( FI +2#*F2+2#*F3+F4 ) 

N=1 

N I S-N I S+ 1 

C INCREASE STEP SIZE IF INTEGRATION SMOOTH ENOUGH 

I F ( DS TEST • GT • • 005 • OR • ABS ( S— S I ) *LE • • 1 E— 05« OR« LLL# GT • 1 )GO TO 55 
DS = 1 • 9* DS 
GO TO 55 

C DELO TOO LARGE - HALVING MODE IF LOWER BOUND IS KNOWN 
530 CELT AU= DELO 

DELO = • 5* (DELTAU+DELTAL ) 

I F ( ABS ( DEL O-DEL TAL ) • LE • • I E-05« AND*LLLL • EQ • 1 ) GO TO 5300 
C SET LLLL TRIGGER - UPPER AND LOWER BOUNDS FOR DELO FOUND - START HALVING MODE 
LLLL=2 
GO TO 2 

C SEARCH AGAIN FOR LOWER BOUND FOR DELO 
5300 DELO=.5*DELO 
DELTAL=DELO 
GO TO 2 

53 I F ( ABS ( UO— • 999*A0ST ) * LE • • 0000 1 ) GO TO 60 
IF ( UO— • 999* AOS T >60*60*54 

C HALVE STEP SIZE TO ACHIEVE VELOCITY BRACKET 

54 DS = DS/2 • 

UO=HUO 
LLL = 2 

55 IF (S*EQ*0*0) GO TO 200 
RS = S 

RUO=UO 

RBETA=BETA 

RDELT A = DELT A 

RRO=RO 

RXO=XO 

RF = F 

GO TO 30 

C VELOCITY WITHIN BRACKET - - PREPARE TO EXTRAPOLATE 
60 S=S+DS 

BETA=BETA+1 ./6.* (Bl +2 • *B2+2 • *B3+B4 ) 
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DELTA = DELTA+1 ./6.* (Cl+2.*C2+2.*C3+C4 ) 

RO=RO+ 1 ./6** (D1 +2.*D2+2.*D3+D4 ) 

X0=X0+1 ./6.* (El +2 .*E2+2.*E3+E4 > 

F = F+ 1 ./6.* (FI +2.*F2+2.*F3+F4 ) 

N'l S = NIS+1 

CALL BSR (LiLLiSI iS ) 

CALL SHOCK ( S * UO * BET A * DELT A « RO ♦ XO * F ) 

IF (L.GT. 1 .AND.PZ1PB.LT. .1 >GO TO 530 
IF (R*DU0DS.LT. • 1 ) GO TO 530 

IF (LN.EQ.O ) WRI TE ( 6 ♦ 99 ) S .RBE TAD * THE T AD «XOi RO * UO * PO i RHOO »T0» DELTA ♦ X 1 
1 iRl *V1 *U1 ,P1 ♦RH01 *DS*DBDS«DUODS*POBAR,SBAR 
IF (SIN (R8ETA )**2.LE. .00002 • OR. ABS ( COS ( A LAMB > > #LE. .00002 .OR.FAC.LT. 
1 O •O.OR . ABS ( RHOO ) .LE. .0000 1 .OR. ABS ( RRO ) .LE • .00001 .OR. ABS (DENOM ) .LE. 
2.00001 .OR # ABS (DEN0M1 > .LE. .00001 )G0 TO 71 
C EXTRAPOLATION TO SONIC POINT 
SST = S+ ( AOST-UO ) / (2.*DU0DS ) 

UOST = AOST 

BETAST=BETA+DBDS* (SST-S ) 

DELST=DELTA+DODS*(SST-S ) 

ROST = RO + DROOS* (SST-S ) 

XOST=XO+DXODS* (SST-S) 

FST=F+DFDS*(SST-S ) 

DS=SST— S 
N I S=N I S + 1 

CALL BSR(L.LL.SI •SST) 

CALL SHOCK ( SST ♦ UOST .BET AST . DEL ST ♦ ROST * XOST ♦ FST ) 

IF (LN.EQ.O ) WRI TE (6*99 > SST * RBET AD ♦ THET AD ♦ XOST . ROST * UOST * PO ♦ RHOO * TO • 
lDELSTiXl «R1 *V1 «U1 «P1 * RHO 1 ♦ DS % OBDS . DUODS « POBAR. SBAR 
IF (LN.EQ.O ) WRI TE (6. 98 ) 

IF (LN.EQ.O )G0 TO 601 
IF (L.EQ. 1 >G0 TO 600 
C TEST FOR SONIC °0 1 NT LOCATION 

IF ( ABS CROST/R-1 • ) .LE. .0001 ) GO TO 600 
IF (ROST.GT.R )GO TO 530 

IF (ROST *LT .R * AND . LLLL • EQ* 1 ) GO TO 602 
IF (R0ST.LT.R.AND.LLLL.EQ.2 ) GO TO 603 
C CORRECT DELO FOUND - - SET LENGTH SCALE AND RERUN FOR PRINTING 

600 LN=0 

CD=4.*FST/R0ST**2 

I F ( L • EQ • 3 ) GO TO 604 

I F ( L . EQ . 1 ) DELO = DELO/SST 

IFCL.EQ.l >DSC=DSC/SST 

SMAX=SST 

GO TO 2 

601 WR I TE ( 6 ♦ 1 00 ) 
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GO TO 71 

C DELO TOO SMALL - - SEARCH AGAIN FOR UPPER BOUND FOR DELO 

602 DELTAL=DELO 
DEL0=2.*DEL0 
GO TO 2 

C DELO TOO SMALL - HALVING MODE 

603 DELT AL = DELO 

DELO= • 5* < DELTAL + DELTAU ) 

GO TO 2 

C RESCALING SPHERE CONE TO A UNIT BASE RADIUS 

604 L=4 

DELO = DELO /R 
DSC=DSC/R 
S I =SI/R 
SMAX=SST/R 
GO TO 2 

70 WR I TE ( 6 « 94 ) ERR ( 1 ) * DUODS 

71 CONTINUE 
GO TO 1 
END 


FUNCTION COT (A) 
COT = COS ( A )/SIN( A) 
RETURN 
END 


FUNCTION TAN<A) 
TAN=SIN(A)/COS( A) 
RETURN 
END 


SUBROUTINE BSR ( L ♦ LL * S I t S ) 

COMMON ER ( 12 ) ,BL< ( 12 ) *6CK< 12 ) • ERR (12) • THETA 
1 • C ON 1 < C0N2 *C0N3 «CON4 « CON5 * CON6 * CON7 * C90 *C91 
21*STlfZl«Hl*Rl*Xl * CPH I * RHOO tPO ♦ AOSQ » STO «PH I 
3DXODS • TO • AOST »BETA»FACi DENOM < DENOM 1 ,DFDSiCD 
4SBAR,LNiDELTAL t DELO • DELT AU # N I R * N I S * L « PZ 1 PB 
GO TO( 1 »2,3<4,5)iL 
C DISK 

1 THET A= 1 .5707963 
AK = 0.0 
GO TO 9 


♦ THET AD « AK ♦ RBET AD . AL AMB 

♦ NNN.U1 .VI «P 1 ♦ GAMMA « RHO 
* DDDS * DBDS * DUODS .DRODS ♦ 
« R i PU 1 PB . POBAR i SMAX * 
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C CONVEX SPHERICAL CAP PORTION - - LENGTH SCALE IS SPHERE RADIUS 

2 THETA=1 .5707963-S 
AK = 1*0 

GO TO 9 

C SPHERE-CONE CHECK FOR JUNCTION 

3 IFC (SI— S).GT. .00001 )G0 TO 2 
IF(ABS(S-SI > .LE. .00001 .AND. LL.EQ. 1 ) GO TO 6 

C SPHERE-CONE BEYOND JUNCTION 
AK = 0 • 0 
GO TO 9 

C RESCALED SPHERE-CONE CHECK FOR JUNCTION 

4 IF ( < S I — S ) .GT ♦•00001 ) GO TO 7 

IF (ABS (S-SI ) *LE. *00001 •AND.LL.EQ* 1 )G0 TO 8 
C RESCALED SPHERE-CONE BEYOND JUNCTION 
AK=0.0 
GO TO 9 

C CONCAVE SPHERICAL CAP 

5 THETA=1 .5707963+S 
AK = — 1 .0 

GO TO 9 

C SPHERE-CONE JUNCTION REACHED 

6 LL = 2 

GO TO 2 

C RESCALED SPHERE-CONE ON SPHER I CAL PORTION - LENGTH SCALE IS CONE BASE RADIUS 

7 THETA=1 .5707963— R*S 
AK =R 

GO TO 9 

C RESCALED SPHERE-CONE JUNCTION REACHED 

8 LL =2 

GO TO 7 

9 THET AD = THET A *57 . 295780 
RETURN 

END 


SUBROUTINE SHOCK (RSiRUO ♦ RBETA * RDELT A . RRO .RXO.RF) 

COMMON ER < 12 ) . BLK (12) . BCK (12) ♦ ERR < 1 2 ) . THETA « THET AD . AK . RBET AD • ALA MB 
1 .CONI .CON2.CON3.CON4.CON5.CON6.CON7.C90.C9I .NNN.Ul .VI »P1 . GAMMA. RHO 
2I.ST1.Z1.H1.R1.X1 « CPH I . RHOO . PO . AOSG . STO , PH I . ODDS . DBDS . DUODS . DROD S • 
3DXODS ♦ TO . AOST . BET A, FAC . DENOM . DENOM 1 , DFDS ♦ CD . R . PU 1 PB . POB AR . SMAX ♦ 
^SBAR.LN. DELTAL . DELO . DELTAU .NIR.NIS.L.PZ1PB 
THET AD= THET A *57. 295780 
RBET AD=RBET A *57. 295780 
C90= 1 .+ak*rdelta 
C91 =2.+AK*RDELTA 
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C SHOCK WAVE FUNCTIONS 
alamb=rbeta-theta 
B=CON2*S I N CRBETA > **2 
A = • 5*CON 1 * ( B— 1 * J/CON2 

I F ( NNN • EQ • 0 ) WR I TE ( 6 « 94 ) BLK ( 8 ) • RS * RUO * RBET A • RDELT A « RRO «RXO* R8ET AO ♦ T 
1HETA*AK» C90*C91 »ALAMB»B*A 

IF(SIN(RBETA)**2.LE..00001.0R.ABS(COS(ALAMB) ).LE..OOOOl )GO TO 1 
U1 = ( 1 .-A )*COS ( THETA >+A*COT < RBET A ) *S I N < THETA ) 

V 1 =— ( 1 • — A } *S I N ( THETA )+A*COT (RBET A )*COS (THETA ) 

P l = ( 2 . *GAMMA*B-CON3 ) / ( GAMMA *CON4*CON2 ) 

RHOl =C0N4*6/ (2 • +CON3*8 ) 

C ISENTROPE CONSTANT FOR NORMAL SHOCK 

I F ( ABS<RS ) *EQ.O .0 )CPH I = P 1 /RHOl **GAMM A 
ST 1 =RHO!*Ul 
Z1 ~ V 1 *ST 1 
H 1 =RHO 1 * V 1 

R 1 =RRO+RDELTA*COS< THETA ) 

X 1 =RXO-RDELTA*S IN ( THETA > 

I F ( NNN • EQ • 0 ) WR ITE(6*94) BLK ( 9 ) * U 1 *Vl «Pl *RH01 « ST1 • Z1 * H 1 • R l « X 1 
FAC=(CON5*< 1 •~RUO**2 )*fCON6 )/CPH I 
IF(FAC.LT.O.O)GO TO 2 

C I SENTROP I C SURFACE PRESSURE AND DENSITY 
RHOO-F AC**C0N7 
PO=CPHI *RHOO** GAMMA 
I F ( RS • EQ • 0*0 )PMAX=PO 
POBAR=PO/PMAX 
I F < LN • EQ • 1 )SMAX=1.0 
SBAR=RS/SMAX 

IF (ABS (RHOO) .LE. .00001 ) GO TO 3 

AOSQ=GAMMA*PO/RHOO 

TO=CON2*AOSQ 

STO=RHOO*RUO 

C SHOCK PARTIAL DERI VS WRT BETA 

PU1 PB = -CONl *COS (RBET A ) * S I N ( ALAMB ) -A*S I N ( THET A ) /SIN (RBET A )**2 
PVlPB = CONl*CDS (RBET A >*COS( ALAMB ) -A*C0S (THETA )/SIN( RBET A ) **2 
PPlPB = CONl *S IN (RBET A ) *COS ( RBET A ) 

PROl PB=RHO! **2*4.*COT (RBET A )/ (CON4*8 ) 

PZ 1 PB=H 1 *PU 1 PB+ST 1 *PV 1 PB+U 1 *V 1 *PR01 PB 
IF (L.GT. 1 .AND.PZ1PB.LT. .1 ) GO TO 7 
DEN0M5=RH0I *CPHI 

I F ( ABS (RHOl ) .LE. .0000 1 .OR. ABS (PI ) .LE. .00001 • OR . ABS ( DENOM5 ) #LE. . 000 
101 ) GO TO 8 
D I SC 1 2 P 1 /DENOM5 
IF (DISCI .LT.0.0 )GO TO 9 
TAU1 -D I SC 1 **C0N7 
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TAUO=RHOO 

PTAUPB=TAUl /CON3* (PPl PB/P1 -PROl PB/RHOl ) 

PU 1 T 1 B=U1 *PTAUPB+TAU1 *PU1 PB 
IF ( RS • EQ • 0 • 0 ) GO TO 7 
IF(ABS(RR0).LE#1*E-13)G0 TO 4 
DERIVATIVES WRT S 

DDDS=C90*TAN ( ALAMB ) 

PHI =RDELTA*S IN (THETA )/RRO 

IF (NNN.EQ.O > WR I TE (6i 94 > BLK ( 10 ) • FAC * RHOO • PO • AOSQ * TO * STO * PH I 
C92=Rl/RRO 

DEN0M=RDELTA*C92*PZ 1 PB 

IF (NNN.EQ.O ) WRI TE (6 .94 ) BLK (11) , PU 1 PB . P V 1 PB . PP 1 PB . PRO 1 PB • PZ 1 PB . DDDS 
1 ♦ C92 . DENOM 

I F ( ABS ( DENOM ) • LE • • 0000 1 )GO TO 5 

FAC1=C90*(TAN(ALAMB)-PHI ) *Z 1 -C9 1 *C92* V 1 *H 1 + AK*RDELT A*RUO*STO+ < C90+ 
1 C92 ) * (PO-P1 ) 

DB0S=FAC1 /DENOM 

DENOM 1 =RDELTA*TAUO* ( AOSQ-RUO**2 ) 

IF (NNN.EQ.O )WR I TE < 6 . 94 ) BLK ( 12) . FAC 1 . DBDS • DENOM 1 iRHOl .PI .RHOO .01 SCI 
1 . T AU 1 tTAUO .PTAUPB.PUIT1 B 
DRODS=S IN (THETA ) 

DXODS=COS (THETA ) 

DFDS= < PO— CONS ) *RRO*S I N ( THETA ) 

I F ( ABS ( DENOM 1 ) • LE **00001 > GO TO 7 

FAC2=C92* (RDELTA*PU1T1B*DBDS+C91*TAU1*V1 )+ (RUO*TAUO— Ul *T AU 1 )*DDDS + 
1PHI*(RUO*TAUO+C90*U1*TAU1 ) 

DUODS=-AOSQ*FAC2/ABS(DENOMl ) 

GO TO 7 

1 WRITE (6 *94 )ER(6) . RBET A « ALAMB 
GO TO 7 

94 FORMAT ( 1X1A6/(8E16.8> ) 

2 WRITE(6. 94 ) ER ( 7 ) . FAC 
GO TO 7 

3 WR I TE < 6 . 94 )ER (8 ) . FAC . RHOO • PO 
GO TO 7 

4 WR I TE ( 6 . 94 )ER(9).FAC. RHOO . PO . AOSQ . TO • STO tRRO 
GO TO 7 

5 WR ITE (6.94 )ER ( 1 0 ) . DENOM . PZ 1 PB . DELTAL . DELO • DELTAU 
WRI TE (6*91 )N IRANIS 

91 FORMAT(IOIS) 

GO TO 7 

8 WR I TE ( 6 . 94 )BCK(6) . RHO 1 .PI .DENOM 5 
GO TO 7 

9 WRITE<6*94 )BCK (7) •DISCI 

7 IF (NNN.EQ.O ) WR I TE ( 6 ♦ 94 ) BCK ( 1 ) * FAC2 . DUODS * DRODS * DXODS • DFDS . CD 
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RETURN 

END 


Input Cards for Sample Cases Described in Appendix B 


ERl 

ER2 

ER3 

ER4 

ER5 

ER6 

ER7 

ER8 

ER9 

ER10 

ERl 1 

ERl 2 

BLK1 

BLK2 

BLK3 

BLK4 

BLK5 

BLK6 

BLK7 

BLK8 

BLK9 

blkio 

BLK1 1 

BLK12 

BCK1 

BCK2 

BCK3 

BCK4 

BCK5 

BCK6 

BCK7 

BCK8 

BC<9 

BCK 1 0 

BCK 1 1 

BCK 12 

ERR 1 

ERR2 

ERR3 

ERR4 

ERR5 

ERR6 

ERR7 

ERR8 

ERR9 

ERR 10 

ERR1 1 

ERR 1 2 


1 

1 1 

14000000+01 

2 2 

1 4000000+0 1 
14000000+01 
3 2 

1 4000000+0 1 
14000000+01 
5 1 

14000000+01 


46300000+01 

46300000+0 1 
46300000+0 1 

46300000+0 1 
46300000+0 1 

47600000+01 


1 0000000+01 

1 0000000+0 1 
1 0000000+0 1 

60000000+02 

60000000+02 

1 0000000+0 1 


10000000+01 

25882000+00 

50000000+00 

40000000+01 

50000000+02 

25882000+00 
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Figure 15.- Effect of 0 C on pressure distribution for a sphere-cone, y = 1.4; IVU, = 10; R n /Rb = 0.25. 
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Figure 18.- Simulated pressure distribution in plane of symmetry for 60° sphere-corn 
obtained for windward and leeward sides by respectively adding or subtracting angle 
the locations of stagnation and sonic points (see text). 
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Figure 19.- Simulated pressure distribution in plane of symmetry for 60° sphere-cone at 5° and 10° angle of attack, y = 1.4; Ma> = 2.96; R n /Rb = 0*25. Calculated results 
obtained for windward and leeward sides by respectively adding or subtracting angle of attack from nominal (60°). Then s-coordinates are shifted and stretched to correct 
the locations of stagnation and sonic points (see text). 
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Figure 20.- Dependence of sonic singularity location on shock-wave standoff distance, y = 1.4; = 10; 9 C = 52°. 
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Figure 21.- Chart of sonic condition for pointed cone. 




National Aeronautics and Space Administration 
Washington, D. C. 20546 


POSTAGE AND FEES PAID 
NATIONAL AERONAUTICS ANE 
SPACE ADMINISTRATION 


OFFICIAL BUSINESS 


FIRST CLASS MAIL 


02U 001 26 51 3DS 6 
AIK FORCE WEAPONS LABDKATC 
kirtland air force BASE, m- 


Tr KISS MADELINE F. CANOV 

T R O f\ O V / l* ( T i 


POSTMASTER: 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 


“ The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof .” 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and Notes, 
and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



